Final  Performance  Report 
for  Optimization  Problems  in  Multisensor  and 
Multitarget  Tracking 

AFOSR  Grant  Number  FA9550-04- 1-0222 

Aubrey  B.  Poore 
Numerica  Corporation 
4850  Hahns  Peak  Drive,  Suite  200 
Loveland,  Colorado  80538 
Email:  Aubrey.  Poore@Numerica.  US 
Phone  Number:  970-461-2000 

Submitted  to 

Dr.  Donald  Hearn,  Program  Manager 
Optimization  and  Discrete  Mathematics 
Air  Force  Office  of  Scientific  Research  /NL 
875  North  Randolph  Street  Suite  325 
Suite  325,  Room  3112  824A 
Arlington,  VA  22203 
email:  donald.hearn@afosr.af.mil 
ph  (703)  696-1142  fax  (703)  696-8450 

February  25,  2008 


Distribution  Statement:  Government  Purpose  Rights. 

The  Government  may  use.  modify,  reproduce,  release,  perform,  display  or  disclose  these  data  within  the  Government  without  restriction,  and  may  release  or  disclose 
outside  the  Government  and  authorize  persons  to  whom  such  release  or  disclosure  has  been  made  to  use.  modify,  reproduce,  release,  perform,  display  or  disclose 
that  data  for  United  States  Government  purposes,  including  competitive  procurement. 

Agreement  Number:  FA9550-04-1-0222 
Recipient:  Dr  Donald  Hearn,  Program  Manager 
Recipient’s  Address: 

Optimization  and  Discrete  Mathematics 
Air  Force  Office  of  Scientific  Research  /NL 
875  North  Randolph  Street  Suite  325 
Suite  325,  Room  3112,  824A 
Arlington,  VA  22203 


1 


200805271 61 


I 


AFRL-SR-AR-TR-08-0276 


REPORT  DOCUMENTATION  PAGE 

Public  reporting  burden  fa  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  ir 
data  needed  and  completing  and  renewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  othei 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Rqoorts  (0704 -< 

4302  Respondent  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  ... 
valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. _ 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

25-02-2008 _  Final  Performance  Report _ 

4.  TITLE  AND  SUBTITLE 


Optimization  Problems  in  Multisensor  and  Multitarget  Tracking 


6.  AUTHOR(S) 

Aubrey  B  Poore 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 


wnii  a  collection  of  information  if  it  does  not  display  a  currently 


3.  DATES  COVERED  (From  -  To) 

9  Jan  04-30  Nov  07 

5a.  CONTRACT  NUMBER 

FA9550-04-1  -0222 

5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


Numerica  Corporation 

PO  Box  271246 

Fort  Collins,  CO  80527-1246 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

'^7‘Ta7  OurloLph  St  3 2 s' 

tlBUTION  /  AVAIL -  ~ 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT:  Government  Purpose  Rights 

The  Government  may  use,  modify,  reproduce,  release,  perform,  display  or  disclose  these  data  within  the  Government  without 
restriction,  and  may  release  or  disclose  outside  the  Government  and  authorize  persons  to  whom  such  release  or  disclosure 
has  been  made  to  use,  modify,  reproduce,  release,  perform,  display  or  disclose  that  data  for  United  States  Government 
purposes,  including  competitive  procurement.  ^  I)  l  ( 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

The  objective  of  this  research  program  is  to  develop  optimization  algorithms  that  solve  key  problems  in  multiple  target  tracking  and  sensor 
data  fusion.  The  central  problem  in  multiple  target  tracking  is  the  data  association  problem  of  partitioning  sensor  reports  into  tracks  and 
false  alarms.  New  classes  of  data  association  problems  have  been  formulated  and  initial  algorithms  developed  to  address  cluster  tracking, 
merged  measurements,  and  even  sensor  resource  management  in  the  form  of  "group-assignments.”  In  a  different  direction,  an  efficient  k- 
best  algorithm  has  been  developed  to  approximate  the  uncertainty  in  data  association,  which  is  critical  for  discrimination  or  combat 
identification.  Statistical  Monte  Carlo  methods  are  also  applicable  and  are  still  under  investigation.  Bias  estimation  algorithms  using  known 
data  association  such  as  truth  objects  and  targets  of  opportunity  have  been  developed.  Bias  estimation  in  which  data  association  is 
unknown  is  difficult  due  to  the  nonconvex  and  mixed  integer  nature  of  the  mathematical  formulation.  Exact  and  approximate  algorithms  have 
been  developed  and  successfully  applied  to  system  tracking.  As  a  prerequisite  to  the  development  of  multiple  target  tracking  approaches  to 
space  surveillance,  consistent  measures  of  uncertainty  for  initial  orbit  determination  and  the  propagation  of  the  uncertainty  over  time  have 
been  developed. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Aubrey  B  Poore 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

None 

31 

including 
this  page 

19b.  TELEPHONE  NUMBER 

(include  area  code) 

(  970)  461-2000 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


1  OBJECTIVES 


1  OBJECTIVES 

The  objectives  of  this  research  program  are  to  develop  optimization  algorithms  that  solve 
key  problems  in  multiple  target  tracking  and  sensor  data  fusion.  While  these  optimization 
problems  are  motivated  by  multiple  target  tracking  applications,  the  class  of  problems 
and  corresponding  algorithms  are  abstracted  to  a  more  mathematical  setting  so  that  the 
resulting  research  developments  have  a  wider  applicability. 

The  first  class  of  these,  the  data  association  problem,  is  the  central  problem  in  mul¬ 
tiple  target  tracking.  This  is  manifest  in  partitioning  measurements  into  tracks  and 
false  alarms,  and  also  in  associating  tracks  from  multiple  sources.  For  measurement- 
to-measurement  or  measurement-to-track  association,  the  multidimensional  assignment 
problem  governs  the  association  process.  Lagrangian  relaxation  are  combined  with  an 
exact  method  such  as  A*-search  or  branch  and  bound,  which  if  run  to  completion,  would 
produce  an  optimal  answer.  (Optimal  solutions  are  required  to  run  off-line  in  simulations 
to  verify  that  the  approximate  algorithm  is  solving  the  association  problem  to  within  the 
noise  level  in  the  problem.)  While  the  basic  algorithms  for  constructing  a  solution  are  now 
mature  for  single  assignments  and  work  exceptionally  well,  new  problems  and  demands 
require  research  into  new  formulations  and  algorithms  supported  by  mathematical  theory. 

The  first  of  these  arise  from  the  need  to  perform  group-formation  tracking  for  ground 
targets,  group-cluster  tracking  and  pixel-cluster  tracking  for  IR  sensors,  and  to  process 
merged  measurements  for  in  (narrowband)  radar.  Thus,  an  objective  is  to  formulate  a 
class  of  new  problems  within  the  framework  of  the  multidimensional  assignment  problem, 
namely  the  “group-assignment”  problems,  which  have  strong  similarities  to  multiple  di¬ 
mensional  combinatorial  auctions.  These  problems  reduce  to  the  traditional  one-to-one 
multidimensional  assignment  when  the  groups  do  not  overlap.  In  addition,  these  same 
group-assignment  problems  also  appear  to  have  much  broader  application  to  new  prob¬ 
lems  arising  in  network  management,  procurement,  or  resource  scheduling  under  the  name 
“combinatorial  auctions.” 

In  addition  to  this  new  class  of  combinatorial  optimization  problems,  an  important 
new  area  for  multiple  target  tracking  is  that  of  supplying  some  measure  of  uncertainty  in 
the  data  association  process.  Here  is  the  reason.  Most  tracking  systems  already  supply 
a  covariance  matrix  to  represent  the  uncertainty  in  the  track  state  due  to  measurement 
noise  and  mismodeling  of  the  dynamics  of  the  object.  This  uncertainty  does  not  represent 
the  data  association  uncertainty,  i.e.,  the  measurement-to-track  or  track-to-track  assign¬ 
ment.  Almost  all  target  identification,  combat  identification,  or  discrimination  algorithms 
assume  perfect  data  association  and  are  not  robust  to  misassociation;  however,  tracking 
systems  cannot  produce  perfect  data  association  given  the  stochastic  nature  of  the  prob¬ 
lem.  Thus,  in  order  to  properly  process  kinematic  and  feature  data  in  the  identification, 
combat  identification,  and  discrimination  process,  an  assessment  of  the  association  pro¬ 
cess  is  required.  This  capability  will  be  required  by  all  future  tracking  systems.  While 
we  have  developed  an  efficient  k-best  algorithm  to  approximate  this  uncertainty,  Markov 
Chain  Monte  Carlo  methods  are  being  investigated  as  well. 

A  third  component  of  this  research  program  is  the  treatment  of  sensor  biases  and 


2 


Approved  for  public  release; 
distribution  is  unlimited. 


1  OBJECTIVES 


navigation  errors,  which  is  a  prerequisite  to  multi-sensor  tracking  and  fusion.  These 
errors  contain  a  deterministic  component  (e.g.,  a  mean)  as  well  as  a  purely  stochastic 
component.  One  can  estimate  a  deterministic  component  using  a  nonlinear  maximum  a 
posteriori  (MAP)  estimation  formulation,  which  is  posed  as  a  weighted  nonlinear  least 
squares  problem  when  the  errors  are  assumed  to  be  Gaussian.  The  least  squares  problem 
requires  training  data  that  can  be  either  from  known  data  association  such  as  targets 
of  opportunity  in  which  the  target’s  dynamic  state  must  also  be  estimated  or  certain  or 
uncertain  truth  data  which  does  not  need  to  be  estimated. 

A  fourth  component  is  the  problem  of  joint  bias  estimation  and  data  association  using 
the  target  themselves  in  which  data  association  is  unknown.  This  leads  to  a  class  of 
nonconvex  mixed  integer  nonlinear  programming  problems.  Once  considered  too  difficult 
for  practical  applications,  this  problem  is  taking  center  stage  due  to  the  need  to  correct 
biases  in  track  states  transmitted  over  a  network  where  access  to  the  fine  details  of  the 
sensor  measurement  process  is  not  available.  This  problem  is  central  to  processing  UAV 
track  states,  space  surveillance,  system  level  tracking,  and  missile  defense. 

A  fifth  component  of  this  work  has  been  the  initial  investigations  into  initial  orbit 
determination,  the  uncertainty  associated  with  the  estimates,  and  the  propagation  of  the 
uncertainty  over  time.  A  separate  report  is  attached  to  this  one  and  is  summarized  below. 

While  the  above  four  classes  of  problems  remain  a  focus,  Numerica  continues  the  work 
on  network-centric  tracking  that  was  initially  started  with  AFOSR  funding  and  has  now 
transitioned  to  a  Phase  II  SBIR  at  AFRL/SNAT  and  a  Phase  II  SBIR  with  Department 
of  the  Army  with  a  transition  path  to  the  SIAP  Joint  Program  Office. 

1.1  STATUS  OF  EFFORT 

The  status  of  the  current  research  program  is  summarized  in  this  section.  The  individual 
topics  are  addressed  in  the  following  subsections.  Of  these,  the  most  significant  part  of 
the  effort  in  the  last  year  has  been  the  development  of  a  general  approach  to  the  joint 
bias  estimation  and  data  association  discussed  below. 

1.1.1  The  Group- Assignment  Problem 

Over  the  last  fifteen  years,  the  multidimensional  assignment  problem  in  which  one  assigns 
each  measurement  or  track  at  most  once  has  generally  been  accepted  as  the  correct  for¬ 
mulation  of  the  central  data  association  problem  for  multiple  hypothesis  tracking  (MHT) 
and  is  now  the  industry  standard.  (This  does  not  mean  that  all  systems  have  adopted 
it;  many  US  military  systems  use  very  old  tracking  technologies.)  Improvements  in  the 
existing  algorithms  and  the  development  of  new  algorithms  for  this  problem  will  remain 
a  fundamental  research  problem  for  decades  to  come. 

While  the  current  research  program  has  pursued  improvements  in  these  algorithms,  a 
focus  has  been  to  investigate  new  classes  of  data  association  problems  arising  from  a  need 
to  extend  current  capabilities  to  tracking  large  numbers  of  closely  spaced  objects,  tracking 
in  dense  clutter  on  the  one  hand,  breaking  up  merged  measurements  over  what  a  radar 
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or  IR  sensor  can  normally  resolve.  These  represent  three  different  classes  of  problems: 

a)  one  in  which  measurements  are  merged  or  clustered  to  improve  run-time  performance, 

b)  individual  object  tracking  in  which  each  object  is  tracked,  and  c)  measurements  or 
clusters  are  decomposed  to  better  track  subclusters  or  individual  objects.  In  a)  and  c), 
one  can  use  clustering  methods  to  formulate  the  data  association  problem.  We  now  give 
a  brief  description  of  these  problems. 

The  idea  behind  problem  classes  a)  and  c)  is  to  use  clustering  methods  to  either 
group  measurements  together  or  to  break  up  clusters  or  merged  measurements.  If  the 
correct  number  of  clusters  is  known  a  priori,  then  the  problem  could  be  posed  as  a  normal 
assignment  problem  that  allows  multi-assignment  to  split  or  merge  clusters.  Determining 
the  number  of  clusters  in  cluster  tracking  or  in  clustering  methods  in  general  is  a  key 
problem.  While  there  are  many  suggested  methods  for  determining  the  number  of  clusters 
given  just  one  look  at  the  data,  most  of  these  are  not  robust  in  tracking  where  one  makes 
multiple  looks  at  the  data  (frames  of  data)  as  it  evolves  in  time.  To  take  advantage  of 
this  feature  of  the  tracking  problem,  we  have  proposed  multiple  clustering  of  the  data  on 
each  frame  and  then  deciding  on  the  correct  number  on  a  particular  frame  by  examining 
many  frames  at  once. 

The  corresponding  mathematical  formulation  of  the  problem  can  be  stated  as  follows. 
Let  P  and  Q  denote  two  lists  of  objects  and  let  V  =  {P,},6/  and  Q  =  {Qj}jej  denote 
collections  of  subsets  of  P  and  Q,  respectively.  These  collections  represent  all  the  subclus¬ 
ters  obtained  from  multiple  clusterings  of  the  data.  A  general  formulation  of  the  cluster 
tracking  multi-assignment  problem  allows  the  multi-assignment  between  groups  P*  of  V 
and  Qj  of  Q  directly  and  then  adds  the  set  packing  as  additional  constraints. 


Minimize  L  C%j  Xij  f 

(i,j)€A 

Subject  To  T  %ij  <  mi  (i  G  /), 

JeA(i) 


for  which  ii  ^  i2  and  Pn  n  Pl2  ^  0 
or  ji  ±  j2  and  Qn  n  Qn  ^  0, 


(1) 


G  {0,  l}. 


The  constraint  (1)  is  the  constraint  on  the  (hard)  set  packing  requirement  in  clustering. 
It  essentially  says  that  a  measurement  cannot  be  in  two  different  clusters  in  the  final 
assignment. 

More  extensive  discussions  of  this  formulation  as  well  as  extensions  to  higher  dimen¬ 
sions  are  presented  in  the  work  of  Gadaleta  and  Poore  [3,  13,  4,  12,  11],  The  formulation 
of  the  cluster  assignment  problem  is  particularly  well-suited  to  a  Lagrangian  relaxation 
algorithm  in  that  the  set  packing  constraint  can  be  Lagrangian  relaxed  to  the  base  prob- 
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lem  consisting  of  either  the  usual  one-to-one  assignment  problem  or  the  multi-assignment 
problem.  The  nonsmooth  optimization  of  the  resulting  problem  is  relatively  straightfor¬ 
ward.  The  final  step  in  restoring  the  set  packing  constraint  is  that  which  remains. 

Clustering  Methods.  Clustering  Methods  is  a  very  large  subject  indeed,  and  there 
are  so  very  many  clustering  methods  [15].  Since  the  error  in  measurements  is  normally 
Gaussian,  we  have  used  a  Gaussian  sum  approximation  and  the  expectation  maximization 
(EM)  algorithm  to  determine  the  means  and  covariances  as  well  as  the  cheaper  K- Means 
algorithm.  A  new  approximation  due  to  Maybeck  and  Williams  [16]  seems  to  be  a  very 
promising  approach  that  can  replace  the  EM  algorithm.  This  has  yet  to  be  investigated 
at  Numerica.  (Note  that  the  track  state  updates  are  clustered  in  the  work  of  Maybeck 
and  Williams  while  the  measurements  themselves  are  clustered  in  this  work.) 

The  Status.  The  following  gives  a  brief  description  of  the  major  tasks  of  this  proposed 
effort.  Though  not  mentioned  explicitly  earlier,  each  of  these  formulations  has  been  tried 
on  small  problems  where  the  optimal  solution  can  be  obtained  by  explicit  enumeration, 
and  the  results  have  been  shown  to  be  very  effective  in  cluster  tracking.  The  number  of 
problems  listed  here  exceeded  our  budget  and  future  work  will  need  to  focus  on  these 
problems  to  further  develop  the  algorithms;  however,  they  are  all  essential  pieces  to  solve 
the  next  generation  of  data  association  problems. 

1.  Two-Dimensional  Multi- Assignment  Solver.  Numerica  has  proposed  the  ex¬ 
tension  of  the  JV  algorithm  to  the  multi-assignment  case  since  it  is  the  fastest  known 
algorithm  for  the  one-to-one  assignment  problem.  An  algorithm  has  been  completed 
that  incorporates  the  enumeration  of  the  K-best  solutions  based  on  Murty’s  method. 

2.  Two-Dimensional  Group-Assignment  Solver  (With  and  Without  Multi- 
Assignments).  The  easiest  problem  in  this  group  is  the  one  without  multi-assignments. 
Numerica  uses  this  formulation  in  tracking  using  UAV  platforms. 

1.1.2  Assignment  Ambiguity  Assessment  and  the  Status. 

The  initial  approach  to  assessing  ambiguity  in  data  association  was  to  compute  the  K- 
Best  solutions  up  to  some  maximum  value  of  K  or  to  some  point  at  which  the  objective 
function  value  drops  below  some  threshold  relative  to  the  optimal  solution.  These  costs 
are  then  used  to  estimate  (optimistically)  the  probabilities  of  association.  While  this  has 
worked  well  for  sparse  to  medium  density  problems,  a  more  robust  algorithm  is  needed 
for  the  medium  to  dense  problems  and  one  that  augments  the  K-Best  solutions.  Work  in 
this  direction  is  proceeding  using  K-Best  and  an  exploration  technique  to  sample  parts  of 
the  solution  space  to  ensure  that  declarations  of  high  probability  of  association  are  indeed 
correct.  Numerica  has  investigated  statistical  Monte  Carlo  Methods,  e.g.,  importance 
sampling  and  Markov  Chain  Monte  Carlo  methods.  The  results  will  be  published  in 
future  work. 
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1.1.3  Bias  Estimation. 

Sensor  bias  and  platform  navigation  errors  represent  systematic  unmodeled  measurement 
and  navigation  errors,  i.e.,  they  cannot  be  averaged  out.  In  single  sensor  environment, 
biases  can  lead  to  optimistic  and  diverging  filters.  In  multi-sensor  environments  bias  errors 
can  lead  to  data  misassociation  (due  to  optimistic  covariances  or  large  residual  biases) 
and  erratic  stereo  multiple  motion  model  performance.  (In  the  worst  case,  N  sensors 
can  produce  N  times  the  true  number  of  tracks.)  There  are  various  classifications  of  the 
biases  such  as  absolute  vs.  relative  and  sensor  biases  vs.  navigation  errors.  In  relative 
bias  estimation,  some  of  the  biases  are  determined  relative  to  others;  the  corresponding 
estimation  techniques  generally  assume  these  to  be  known  a  priori. 

A  key  requirement  for  bias  estimation  is  the  need  for  training  data.  There  are  many 
proposals  in  the  literature  including  targets  of  opportunity  in  which  data  association  is 
known  by  the  various  sensors  and  truth  objects,  e.g.,  precision  reports  based  on  GPS.  In 
general,  the  state-of-the-art  is  to  use  targets  with  known  data  association  spread  over  the 
field  of  regard  to  increase  the  number  of  biases  that  are  observable  from  the  data.  Such 
data  may  not  always  be  present  or  is  very  difficult  to  obtain.  In  the  quest  for  a  rich  set 
of  training  data,  one  could  turn  to  the  targets  themselves,  but  here  the  problem  is  that 
the  data  association  is  unknown.  This  has  been  considered  to  be  too  difficult  a  problem 
in  that  it  results  in  a  nonconvex  mixed  integer  nonlinear  programming  problem.  Thus, 
one  can  delineate  the  bias  estimation  problems  into  two  classes:  a)  those  with  known 
data  association  using  either  truth  objects  or  objects  of  opportunity  and  b)  those  with 
unknown  data  association.  We  now  address  both  of  these  problems. 

The  methods  being  investigated  herein  use  numerical  linear  algebra  techniques  in  the 
estimation  procedure  (e.g.,  singular  value  decomposition)  to  determine  the  relative  biases 
a  posteriori  as  part  of  the  solution. 

Bias  Estimation  Using  Known  Data  Association.  In  funding  from  both  AFOSR, 
AFRL/SNAT,  and  Numerica  internal  funding,  Numerica  has  undertaken  systematic  treat¬ 
ment  of  biases  based  on  “batch  ML/MAP  estimation  with  process  noise,”  which  in  many 
cases  reduces  to  a  nonlinear  least  squares  problem.  (The  correct  data  association  is  as¬ 
sumed.)  The  singular  value  decomposition  is  used  to  determine  the  bias  roll-ups  and  the 
observable  biases.  The  initial  results  [10]  and  later  more  careful  analysis  of  Herman  and 
Poore  [5]  using  truth  objects  and  the  work  of  Kragel,  Herman,  Danford  and  Poore  [7] 
lead  to  performance  of  the  multisensor  tracking  comparable  to  performance  without  bi¬ 
ases.  These  works  are  based  on  a  nonlinear  least  squares  formulation  and  the  the  singular 
value  decomposition  to  determine  the  observable  and  unobservable  biases. 

Here  is  a  brief  technical  description  of  the  problem.  In  the  absence  of  bias,  a  least- 
squares  estimate  for  the  target  state  based  on  noisy  measurements  is  obtained  by  mini- 
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1  OBJECTIVES 


mizing  the  objective  function  (see  Jazwinski[6],  Sage  and  Melsa[14]  or  Poore[9]  et  al.) 

i  Nt  ( 

2E  E  H4S)-4')(4<))ll^(.))-i  (2) 

»=1  \(s,k)€Mi 

+  E  114°  “  A-i(4-i)IIq->  +  Il4°  - 

fc= 2 


where  Tt  is  the  number  of  (unique)  time  updates  for  target  i,  Fk  is  the  Jacobian  of  /*, 
P0  is  the  covariance  of  the  state  prior,  and  s  and  k  denote  the  sensor  index  and  time 
index,  respectively.  The  first  term  in  (2),  the  sensor  component ,  embodies  the  influence 
of  the  observations  on  the  estimated  target  track,  while  the  second  term,  the  process  noise 
component ,  enforces  model  fidelity.  The  third  term,  the  state  prior  component,  incorpo¬ 
rates  information  based  on  previous  estimates.  Note  that  the  Mahalanobis  distance  is 
used  in  (2)  rather  than  the  Euclidian  distance.  This  strongly  influences  the  result  of  the 
estimation  procedure  since  the  relative  influence  of  each  of  the  components  is  determined 
by  the  relationship  between  the  matrices  and  Pq. 

Adding  bias  estimation  to  (2)  is  now  straightforward.  We  replace  the  measurement 
model  in  the  sensor  component  with  the  modified  model  of  (3)  and  append  bias  process 
noise  and  bias  prior  components  similar  to  the  corresponding  state  components.  This 
yields  the  objective  function 


*<x)-jE  E  iK(^.r,',)-4‘,(4',,C(‘l)iiL,)-,  (3) 

»=1  \(s,k)€Mi 

+  E  11**°  -  A-i(4-i)IIq-‘  +  Iki0  -  /i(*o0)ll (fi/vf+QO-*  )  (4) 

k= 2  / 

+  5  E  (E  U4*’  -  +  Ilf’S*’  -  /S(f’!,*,)ii(ffPi,ft,r+0i,-.')  •  (5) 

8=1  \k= 2  / 


where  Ns  denotes  the  number  of  sensors,  Ts  is  the  number  of  (unique)  update  times  from 
sensor  s  and 

X  =  {x£\  &f)(s),  b[M){3) } ,  (s,k)eMi,i  =  1, . . . ,  Nt,  (6) 

is  the  vector  of  unknown  states  and  biases  for  all  targets  and  sensors  in  the  current 
estimation  window. 

Once  sensor  biases  and  navigation  errors  are  estimated  and  the  biases  and  errors 
corrected,  one  is  still  left  with  the  stochastic  part  of  the  biases  called  residual  biases. 
We  have  dealt  with  these  via  an  extension  of  the  Schmidt-Kalman  filter  to  the  nonlinear 
problem  and  multiple  interacting  filter  models  [8].  The  combination  extensively  improves 
tracking  performance. 


Status.  This  work  is  now  complete. 
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Bias  Estimation  Using  Targets  with  Unknown  Data  Association.  The  simplest 
of  these  problems  occurs  in  trying  to  match  track  states  from  two  different  sources.  This 
is  a  key  problem  in  associating  and  fusing  track  states  from  UAVs  or  in  space  surveillance 
where  the  sensor  reports  are  track  states  rather  than  measurements.  The  problem  can  be 
posed  as  the  problem  of  determining  a  displacement  d  and  an  assignment  x  that  optimize 
the  the  mixed  integer  nonlinear  programming  problem 


Minimize (l,d)  Cr(d)  +  Y,  Cij(d)xij 

(*J)eA 

Subject  To:  ^  Xij  <1  (i  = 

j£A(i) 

Y  Xii  -  1  0’  =  !>  •••»«). 

Xij  £  {0, 1} 


(7) 


(8) 


where  A  is  the  set  of  feasible  arcs,  A(i)  =  {j  |  (i,j)  £  A}  and  B(j )  =  {i  |  (i,j)  £  A}  and 
cT{d)  =  ^ dTR~ld 

Cij{d)  =  ^(: ii  +  d-  yjfZ'/ixi  +  d-  yj)  (9) 

+  ^  In  [(3?  det  (27rZjj)]  +  7 jj  for  i  ^  0  and  j  ^  0,  (10) 

Zij  =  Pu~  Sij  —  Sji  +  Qjj  (11) 

where  (JT  denotes  the  target  density,  7 jj  =  hi  ((1  -  Pd)(l  —  P])),  and  P%,  k  £  {1,2}  is 
the  probability  of  sensor  k  observing  a  particular  target. 

In  addition,  there  is  the  multidimensional  assignment  version  of  this  problem  [1], 
(Note  that  for  each  assignment  x,  there  is  a  unique  local  minimum  d.)  The  costs  £Y,(d) 
can  be  general  functions  of  d,  but  the  algorithms  other  than  the  estimation  algorithms 
for  d  given  an  assignment  x  remain  pretty  much  unchanged.  The  first  class  of  these 
algorithms  are  based  on  combined  use  of  heuristics  and  branch-and-bound  algorithms  or 
A*-search.  This  has  been  a  major  achievement  of  the  program  over  the  last  year  with 
preliminary  work  being  reported  by  Danford,  Kragel,  and  Poore  [1,  2]. 


Status.  The  algorithms  [1,  2]  for  the  two  dimensional  problem  are  now  complete,  but 
will  need  to  be  coded  in  an  implementation  language  such  as  C++.  We  have  developed 
four  classes  of  algorithms:  (a)  optimal  algorithm  based  on  branch  and  bound,  (b)  an 
optimal  algorithm  based  on  A*-search,  (c)  a  heuristic  that  is  very  efficient  in  runtime  and 
is  near  optimal,  and  (d)  an  anytime  algorithm.  A  formulation  for  the  multidimensional 
problem  is  complete  and  algorithms  have  been  planned,  but  not  yet  implemented.  This 
work  will  be  submitted  for  additional  publications  in  the  future. 
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2  PERSONNEL  SUPPORTED 


1.1.4  Network-Centric  Multiple  Frame  Assignments. 

An  objective  in  this  research  program  is  the  development  of  distributed  multiple  frame 
data  association  algorithm  that  is  comparable  in  quality  to  that  of  centralized  tracking 
while  managing  communication  loading  and  achieving  a  consistent  or  single  integrated 
air  picture.  This  work  was  initially  funded  by  AFOSR;  however,  due  to  its  importance, 
AFRL/SNAT  and  SIAP  Joint  Program  Office  (SIAP  JPO)  both  funded  Phase  II  SBIRs. 
This  work  is  ongoing  with  strong  potentials  for  transition  to  SIAP  JPO. 

1.1.5  Generation  and  Propagation  of  Track  States  and  Their  Uncertainty  in 
Space  Surveillance. 

The  objective  in  this  research  component  is  the  generation  of  the  initial  orbit,  its  un¬ 
certainty,  and  the  propagation  of  this  uncertainty  over  time.  This  investigation  is  a 
prerequisite  to  the  development  of  an  MHT  tracker  for  space  surveillance. 

Status.  An  attached  report  goes  into  much  more  detail.  This  work  has  just  begun  and 
we  hope  to  continue  in  a  new  AFOSR  contract. 

1.2  ACCOMPLISHMENTS/NEW  FINDINGS 

The  central  problem  in  any  surveillance  system  is  the  data  association  problem  of  par¬ 
titioning  observations  into  tracks  and  false  alarms.  Over  the  last  fifteen  years  and  with 
support  from  AFOSR,  a  new  approach  has  been  developed  based  on  the  use  of  multidi¬ 
mensional  assignment  problem  formulation  and  Lagrangian  relaxation  algorithms.  (This 
approach  is  often  called  multiple  frame  assignments  or  MFA  for  short.)  Four  U.  S.  patents 
have  now  been  issued  for  this  work.  What  is  more,  based  on  this  new  technology,  Lockheed 
Martin  of  Owego,  NY  won  the  best  of  Breed  Tracking  Contest  for  the  next  upgrade  to 
AWACS  held  at  Hanscom  AFB  in  Boston  in  1996,  and  it  has  been  chosen  as  the  tracking 
system  for  the  Navy’s  new  multipurpose  helicopter  under  the  LAMPS  program.  Based 
on  much  of  this  work,  Numerica  has  now  developed  a  system  level  tracker  for  the  Missile 
Defense  Agency  that  should  be  deployed  in  2008. 

The  research  performed  here  is  fundamental  to  future  needs  in  tracking  and  surveil¬ 
lance  over  the  next  fifteen  to  twenty  years. 

2  PERSONNEL  SUPPORTED 

a.  PI:  Aubrey  B.  Poore 

b.  Colleagues  at  Numerica:  Eric  Kuo,  Scott  Danford,  and  Scott  Lundberg. 

c.  Fritz  Obermeyer,  MS  in  Mathematics,  Colorado  State  University,  2005,  Currently 
PhD  student  at  Carnegie  Mellon  University,  Pittsburgh,  PA. 
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3  PUBLICATIONS 


3  PUBLICATIONS 

a.  F.  H.  Obermeyer  and  A.  B.  Poore,  A  Bayesian  Network  Tracking  Database,  in  Oliver  E 

Drummond,  editor,  2004  SPIE  Conference  Proceedings,  Signal  and  Data  Processing 
of  Small  Targets,  SPIE,  August  2004. 

b.  S.  Gadaleta  and  A.  B.  Poore  and  Sean  Roberts  and  B.  J.  Slocumb,  Multiple  Hypothesis 

Clustering  Multiple  Frame  Assignment  Tracking,  In  Oliver  E.  Drummond,  editor, 
Signal  and  Data  Processing  of  Small  Targets,  SPIE  Vol.  5428,  August,  2004,  294- 
307. 

c.  A.  B.  Poore  and  S.  Gadaleta,  The  Group  Assignment  Problem  Arising  in  Multiple 

Target  Tracking,  in  D.  Grundel  and  R.  Murphey  and  P.  Pardalos,  Theory  and 
Algorithms  for  Cooperative  Systems,  World  Scientific,  pages  423-449,  2004. 

d.  Aubrey  B  Poore  and  Sabino  Gadaleta,  The  Group  Assignment  Problem  in  Multiple 

Target  Tracking,  to  appear  in  Proceedings  of  the  4th  International  Conference  on 
Cooperative  Control  and  Optimization,  2005. 

e.  Roman  Novosolev,  Aubrey  B  Poore,  and  Brian  J.  Suchomel,  Network  Centric  MFA 

Tracking  for  Air  Defense:  A  Tribute  to  Oliver  E  Drummond,  2005. 

f.  Sabino  Gadaleta,  Shawn  Herman,  Scott  Miller,  Fritz  Obermeyer,  Benjamin  Slocumb, 

Aubrey  B  Poore,  and  Mark  Levedahl,  Short-term  Ambiguity  Assessment  to  Aug¬ 
ment  Tracking  Data  Association  Information,  Eighth  International  Conference  on 
Information  Fusion,  July,  2005,  Philadelphia,  PA. 

g.  R.  Novoselov  and  S.  Herman  and  S.  Gadeleta  and  A.  Poore,  Mitigating  the  effects  of 

residual  biases  with  Schmidt-Kalman  filtering,  Eighth  International  Conference  on 
Information  Fusion,  July,  2005,  Philadelphia,  PA. 

h.  Aubrey  B  Poore,  Complexity  Reduction  in  MFA/MHT  Tracking,  in  Oliver  E  Drum¬ 

mond,  editor,  2005  SPIE  Conference  on  Small  Targets,  Volume  5913,  August,  2005. 

i.  Allison  Floyd,  Sabino  Gadaleta,  Dan  Macumber,  Aubrey  B  Poore,  Integration  of  IR 

Sensor  Clutter  Rejection  Techniques  with  Pixel  Cluster  Frame  Manipulation,  in 
Oliver  E  Drummond,  editor,  2005  SPIE  Conference  on  Small  Targets,  Volume  5913, 
August,  2005. 

j.  R.  Novoselov  and  S.  Gadeleta  and  A.  Poore,  Network  Centric  MFA  Tracking  Archi¬ 

tecture  Based  on  Soft-Level  Data  Association,  in  Oliver  E  Drummond,  editor,  2005 
SPIE  Conference  on  Small  Targets,  Volume  5913,  August,  2005. 

k.  Daniel  L  Macumber,  Sabino  Gadaleta,  Allison  Floyd,  and  Aubrey  B  Poore,  Hierarchi¬ 

cal  Closely-Spaced  Object  (CSO)  Resolution  for  IR  Sensor  Surveillance,  in  Oliver  E 
Drummond,  editor,  2005  SPIE  Conference  on  Small  Targets,  Volume  5913,  August, 
2005. 
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4  INTERACTIONS/TRANSITIONS 


l.  Shawn  M.  Herman  and  Aubrey  B  Poore,  Nonlinear  Least-squares  Estimation  for  Sensor 

and  Navigation  Biases,  in  Oliver  E.  Drummond,  Proceedings  of  the  2006  SPIE 
Conference  on  Signal  and  Data  Processing  of  Small  Targets.  Vol  6236,  May  2006. 

m.  Benjamin  J  Slocumb  and  Aubrey  B  Poore,  Complexity  reduction  in  MHT /MFA  track¬ 

ing:  Part  II,  hierarchical  implementation  and  simulation  results,  in  Oliver  E.  Drum¬ 
mond,  Ed.,  Proceedings  of  the  2006  SPIE  Conference  on  Signal  and  Data  Processing 
of  Small  Targets.  Vol  6236,  April  2006. 

n.  Aubrey  B  Poore,  Sabino  M  Gadaleta  and  Benjamin  J  Slocumb,  Multiple  Hypothesis 

Correlation  in  TYack-to-Track  Fusion  Management,  in  Gautam  Appa,  Leonidas  Pit- 
soulis,  and  H.  Paul  Williams,  Handbook  on  Modelling  for  Discrete  Optimization, 
Springer  International  Sequence,  Chapter  12,  pp  341-372,  Springer  Verlag,  New 
York,  2006. 

o.  Aubrey  B  Poore  and  Sabino  Gadaleta,  Some  Assignment  Problems  Arising  from  Mul¬ 

tiple  Target  Tracking,  Mathematical  and  Computer  Modelling,  Elsevier,  Volume  43, 
pp  1074-1091,  2006. 

p.  Shawn  M.  Herman  and  Aubrey  B  Poore,  Nonlinear  Least-squares  Estimation  for  Sen¬ 

sor  and  Navigation  Biases,  in  Oliver  E.  Drummond,  Proceedings  of  the  2006  SPIE 
Conference  on  Signal  and  Data  Processing  of  Small  Targets.  Vol  6236,  April  2006. 

q.  Scott  Danford,  Bret  Kragel,  Aubrey  Poore,  Joint  MAP  Bias  Estimation  and  Data 

Association:  Algorithms,  Proceedings  of  SPIE  Conference  on  Small  Targets,  Oliver 
E  Drummond,  Ed.,  Vol  6699,  August,  2007. 

r.  Scott  Danford,  Bret  Kragel,  Aubrey  Poore,  Joint  MAP  Bias  Estimation  and  Data 

Association:  Simulations,  Proceedings  of  SPIE  Conference  on  Small  Targets,  Oliver 
E  Drummond,  Ed.,  Vol  6699,  August,  2007. 

s.  Bret  Kragel,  Shawn  Herman,  Scott  Danford,  and  Aubrey  Poore,  Bias  Estimation  Using 

Targets  of  Opportunity,  Proceedings  of  SPIE  Conference  on  Small  Targets,  Oliver 
E  Drummond,  Ed.,  Vol  6699,  August,  2007. 

t.  Randy  Paffenroth,  Roman  Noveselov,  Marcio  Teixeira,  Stephanie  Chan,  and  Aubrey 

Poore,  Mitigation  of  Biases  Using  the  Schmidt- Kalman  Filter,  Proceedings  of  SPIE 
Conference  on  Small  Targets,  Oliver  E  Drummond,  Ed.,  Vol  6699,  August,  2007. 

4  INTERACTIONS/TRANSITIONS 

a.  Participation/presentations  at  meetings,  conferences,  seminars,  etc. 

i.  Aubrey  B  Poore,  Multiple  Hypothesis  Correlation,  Project  Hercules  Review, 
MITLL  (December  2005  and  November  2006) 
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ii.  Scott  Danford,  Bret  Kragel,  and  Aubrey  B  Poore,  Joint  MAP  Bias  Estimation 

and  Data  Association:  Algorithms  and  Simulations,  SPIE  Conference  on  Signal 
and  Data  Processing  of  Small  Targets,  July  2007. 

iii.  Scott  Danford,  Bret  Kragel,  and  Aubrey  B  Poore,  Joint  MAP  Bias  Estima¬ 

tion  and  Data  Association:  Algorithms,  SPIE  Conference  on  Signal  and  Data 
Processing  of  Small  Targets,  July  2007. 

iv.  Bret  Kragel,  Shawn  Herman,  Scott  Danford,  and  Aubrey  Poore,  Bias  Estimation 

Using  Targets  of  Opportunity,  SPIE  Conference  on  Signal  and  Data  Processing 
of  Small  Targets,  July  2007. 

v.  Randy  Paffenroth,  Roman  Noveselov,  Marcio  Teixeira,  Stephanie  Chan,  and 

Aubrey  Poore,  Mitigation  of  Biases  Using  the  Schmidt-Kalman  Filter,  SPIE 
Conference  on  Signal  and  Data  Processing  of  Small  Targets,  July  2007. 

vi.  Optimization  Problems  in  Multiple  Target  Tracking,  AFOSR  Workshop,  June 

2007,  Washington,  DC. 

v.  Sensor  Resource  Management,  SIAP  Joint  Program  Office,  June  2007. 

b.  Consultative  and  Advisory  Functions 

Dr.  Poore  serves  as  a  Subject  Matter  Expert  on  tracking  and  combat  identification 
for  the  SIAP  Joint  Program  Office  (formerly  JSSEO)  (August,  2003  -  Present). 

c.  Transitions 

i.  Transition:  Numerica 

Numerica  Corporation,  a  Colorado  Corporation,  is  a  small  business  in  Fort 
Collins,  CO  that  is  engaged  in  basic  research,  software  development,  and  en¬ 
gineering  services,  especially  in  surveillance.  Aubrey  B.  Poore  is  President.  In 
January  of  1999,  Numerica  completed  negotiations  with  Colorado  State  Uni¬ 
versity  Research  Foundation  (CSURF)  to  take  out  an  exclusive  license  on  the 
U.  S.  Patents,  software,  and  tracking  technology  developed  by  Aubrey  Poore 
at  Colorado  State  University  (CSU)  (past,  present,  and  future)  for  the  purpose 
of  licensing  the  tracking  technology  to  industry. 

ii.  Transition:  SBIRS  High,  an  Air  Force  Program 

Numerica  participated  in  a  risk  reduction  effort  in  support  of  the  Air  Force 
Program  SBIRS  High.  The  result  of  this  participation  was  a  demonstration 
that  the  tracking  requirements  for  SBIRS  High  could  indeed  be  met. 

iii.  Transition:  MDNTB 

Numerica  transitioned  its  multi- assignment  algorithm  and  software  to  the  Mis¬ 
sile  Defense  National  Team  B  for  use  in  the  ballistic  launch  event  association 
(BLEA)  algorithm  that  will  be  deployed  in  2006. 
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5  NEW  DISCOVERIES,  INVENTIONS,  OR 

PATENTS 


iv.  Transition:  AFRL,  WPAFB 

Numerica  transitioned  its  assignment  solver  and  simulation  environment  for 
urban  surveillance  to  AFRL/SNAT  during  spring  of  2007. 

iv.  Transition:  SBIRS  Low  (Now  STSS) 

Numerica  transitioned  its  MFA  tracking  system  to  Spectrum  Astro  of  Gilbert, 
AZ  (now  General  Dynamics)  and  to  Northrop  Grumman  Corporation  of  El 
Segundo,  CA  for  the  SBIRS  Low  program  (now  STSS). 

v.  Transition:  Project  Hercules,  MDA/BC  and  MDA’s  MDNTB 

In  a  rapid  and  advanced  development  program  funded  by  MDA  for  Project 
Hercules,  Numerica  is  developing  a  state-of-the-art  tracking  system  for  radar 
based  tracking  of  missiles  for  both  national  and  theater  missile  defense.  While 
SBIRS  Low  is  for  satellite  based  IR  sensors,  this  program  focuses  on  radar 
based  systems  located  on  the  ground,  on  the  sea,  and  in  the  air.  As  part  of  this 
effort,  Numerica  has  developed  and  transitioned  the  track  processing  threads 
for  MDA’s  C2BMC  to  MDA’s  National  Team  B.  The  subject  algorithms  are 
being  coded  for  deployment  in  2008. 


5  NEW  DISCOVERIES,  INVENTIONS,  OR 
PATENTS 

Patents  listed  below  to  emanate  from  the  basic  research  supported  by  AFOSR.  Essentially 
an  entirely  new  approach  to  multiple  hypothesis  tracking  (MHT)  has  been  developed. 

a.  NEW  DISCOVERIES/INVENTIONS: 

The  inventions  as  embodied  in  the  U.S.  Patent  below  contains  various  approaches 
to  data  association  problems  based  on  linear  programming  relaxations  and  contains 
a  description  of  the  variable  depth  sliding  window. 

b.  PATENTS  ISSUED  for  AFOSR  SPONSORED  RESEARCH 

Aubrey  B.  Poore,  Jr.,  Method  and  System  for  Tracking  Multiple  Regional  Objects 
by  Multi-Dimensional  Relaxation,  U.S.  Patent  Number  6,404,380  B2,  Issued  11  June 
11  2002.  (Assignee:  Colorado  State  University  Research  Foundation,  Fort  Collins, 
CO) 

Thomas  N.  Barker,  Joseph  A.  Persichetti,  Aubrey  B.  Poore,  Jr.,  and  Nenad  Rijavec, 
Method  and  System  for  Tracking  Multiple  Regional  Objects,  US  Patent  Number 
5406289,  issued  11  April  1995.  (Assignees:  IBM,  Owego,  NY;  Colorado  State  Uni¬ 
versity  Research  Foundation,  Fort  Collins,  CO) 

Aubrey  B.  Poore,  Jr.,  Method  and  System  for  Tracking  Multiple  Regional  Objects 
by  Multi-Dimensional  Relaxation,  CIP,  US  Patent  Number  5537119,  issued  on  16 
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July  1996.  (Assignee:  Colorado  State  University  Research  Foundation,  Fort  Collins, 
CO) 

Aubrey  B.  Poore,  Jr.,  Method  and  System  for  Tracking  Multiple  Regional  Objects 
by  Multi-Diinensional  Relaxation,  CIP,  U.  S.  Patent  Number  5,959,574,  issued  on 
September  28,  1999.  (Assignee:  Colorado  State  University  Research  Foundation, 
Fort  Collins,  CO) 

6  HONORS/ AWARDS 

a.  Name  of  Award:  The  2004  Colorado  State  University  Research  Foundation 
Technology  Transfer  Award  This  award  is  given  yearly  to  one  of  the  faculty  at 
Colorado  State  University  for  achievements  in  transitions  of  technology  to  industry. 
A  complete  description  is  given  at 

http:  /  /  www.csurf.org/ researcher.html. 

Year  Received:  2004 

Honor/Award  Recipient(s):  Aubrey  B.  Poore 

Awarding  Organization:  Colorado  State  University  Research  Foundation 

b.  Name  of  Award:  Colorado  State  University  Alumni  Association  Distin¬ 
guished  Faculty  Award.  (  The  award  is  given  to  a  current  Colorado  State  Uni¬ 
versity  faculty  member,  not  necessarily  a  Colorado  State  University  graduate,  who 
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1  INTRODUCTION 


1  Introduction 

Space  surveillance  is  the  component  of  space  situational  awareness  (SSA)  focused  on  the 
detection  of  space  objects  and  the  use  of  multi-source  data  to  track  and  identify  space 
objects.  Currently,  it  can  take  weeks  to  establish  correct  orbits  on  objects  of  interest, 
due  predominantly  to  the  manual  processes  for  correlation  of  observations  on  an  object. 
This  manual  correlation  process  is  particularly  fragile  when  events  such  as  collisions  or 
the  recent  Chinese  ASAT[2]  test  create  large  numbers  of  new  observable  objects.  In  order 
to  ensure  that  the  space  catalog  can  be  updated  in  a  timely  fashion  when  large  numbers 
of  new  objects  arise,  we  require  modern  multiple  target  tracking  methods. 

Multiple  target  tracking  requires  an  accurate  accounting  of  uncertainty  in  the  measure¬ 
ments,  biases,  and  track  states  for  both  the  initiation  of  tracks  (initial  orbit  determination) 
on  newly  observed  objects,  and  for  the  propagation  of  the  state  estimates  on  those  objects 
over  time.  Currently,  the  space  catalog  does  not  incorporate  state  vector  uncertainty  in 
its  published  Two-Line  Element  (TLE)  sets.  Newer  proposed  data  structures  for  trans¬ 
mitting  information  on  an  object’s  state,  such  as  ISO  22644,  Orbit  Data  Messages  [3], 
may  provide  an  interface  for  communicating  uncertainty  data  on  orbit  states.  While  sig¬ 
nificant  work  has  gone  into  developing  mechanisms  for  high-accuracy  orbit  propagation 
given  an  accurate  initial  orbit  state,  less  focus  has  been  placed  on  the  combined  problems 
of  developing  an  accurate  description  of  the  uncertainty  in  the  initial  orbit  determination, 
and  of  propagating  the  uncertainty  along  with  the  state  vector.  For  example,  for  a  newly 
initiated  orbit,  inaccuracies  in  the  initial  orbit  determination  can  substantially  impact  the 
object  location  days  or  even  hours  later.  The  only  way  to  address  this  inaccuracy  is  to 
develop  an  uncertainty  measure  (e.g.,  a  covariance  matrix)  of  the  new  orbit  that  correctly 
reflects  the  error  in  the  state  vector  estimate.  The  metric  that  we  use  for  describing  the 
correctness  of  the  covariance  matrices  is  “covariance  consistency.”  One  of  the  goals  of  our 
effort  has  been  to  demonstrate  that  consistent  covariances  can  be  achieved  for  the  initial 
orbit  determination  and  that  they  can  be  maintained  when  propagated,  at  least  for  a 
period  of  time  measured  in  days.  The  achievement  of  consistent  covariance  matrices  is  a 
prerequisite  to  the  development  of  an  MHT  tracker  for  space  surveillance. 

Generation  of  consistent  covariances  on  initial  orbit  determination  requires  a  consider- 
analysis  or  Schmidt  filter  to  account  for  residual  biases  or  ill-conditioning  in  the  nonlinear 
transformations  between  state  and  measurement  space.  Propagation  of  uncertainty  over 
time  is  also  challenging,  in  that  propagation  of  the  linearized  uncertainty  in  a  Cartesian 
coordinate  frame  will  not  preserve  consistency.  Instead,  covariance  propagation  requires 
use  of  some  form  of  the  classic  orbital  elements.  The  work  herein  focuses  on  measurements 
from  radar  sensors,  but  the  general  approach  with  appropriate  modifications  can  also  be 
used  for  other  sensors  (such  as  optical  telescopes). 

Section  2  discusses  the  details  of  our  radar  model  and  the  target  motion.  Section  3 
presents  the  algorithmic  approaches  to  achieving  consistent  covariances  for  initial  orbit 
determination  and  for  the  propagation  of  covariances.  Section  4  shows  the  performance 
of  these  methods  in  simulation. 
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2  PROBLEM  STATEMENT 


2  Problem  Statement 

Given  a  sequence  of  radar  measurements  emanating  from  the  same  object,  our  goals  are 
(a)  to  establish  an  initial  orbit,  (b)  to  determine  a  consistent  covariance  matrix  for  that 
object,  and  (c)  to  propagate  the  covariance  or  uncertainty  over  an  extended  period  of 
time  so  that  the  uncertainty  remains  consistent.  Given  the  limited  nature  of  this  effort, 
our  models  are  simplified;  however,  we  believe  the  general  methods  can  be  extended  to 
more  general  force  models  and  different  sensors. 

2.1  Uncertainty  in  Radar  Measurements 

Errors  in  radar  measurements  (range,  azimuth,  and  bearing)  are  significantly  greater 
in  the  angle  components  of  the  measurement  than  in  range.  This  discrepancy  leads  to 
an  error  ellipsoid  that  looks  rather  pancake-shaped  in  Cartesian  coordinates.  Figure  1 
illustrates  this  confidence  region. 


x  10*  % 


°o 


Figure  1:  An  example  uncertainty  region  for  a  measurement.  Each  circle  represents  a 
possible  measurement  the  radar  would  produce  for  a  target  located  at  the  red  star.  From 
this  angle,  the  full  dimensions  of  the  pancake  cannot  be  seen.  Only  one  of  the  two  angle 
errors  is  apparent  here.  This  figure  is  meant  to  demonstrate  the  thin  range  component, 
and  curvature  due  to  nonlinearity. 

In  these  tests,  the  radar  collects  measurements  in  sine-space.  This  space  is  viewed  in 
terms  of  a  cartesian  coordinate  system  aligned  with  the  radar  array  face,  and  originating 
at  the  radar.  The  direction  to  the  target  is  measured  as  the  sine  of  the  angle  between 
the  local  target  vector  and  the  axes  of  the  cartesian  coordinate  system.  The  radar  only 
measures  two  sine  components;  these  two  components  in  combination  with  the  target 
range  suffice  to  compute  the  third  component  if  desired.  If  a  target  is  located  at  position 
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( x ,  y,  z )  in  the  sensor  cartesian  system,  then  the  full  sensor  measurement  is: 

r  =  \/x2  +  y2  +  z2  +  vr 

u  =  -  +  vu 
r 
V 

V  =  -  +  Vv, 

r 

where  uT ,  uu,  and  vv  are  small  independent  errors  that  follow  a  Gaussian  distribution  with 
known  covariance  [5]. 


2.2  Dynamics  Model 


While  some  aspects  of  the  force  models  relevant  to  orbit  dynamics  are  well  characterized 
(such  as  the  gravitational  field),  others  contain  effects  that  may  be  difficult  to  model  ac¬ 
curately.  For  example,  low  Earth  orbit  satellites  may  experience  a  highly  uncertain  drag 
force,  both  because  our  models  of  atmospheric  density  are  inaccurate  and  because  the 
drag  coefficient  of  an  object  will  vary  with  its  orientation  (and  for  many  objects,  particu¬ 
larly  debris  resulting  from  an  event,  it  may  not  be  possible  to  produce  a  drag  coefficient 
estimate).  Additionally,  some  elements  of  the  full  force  model  (such  as  force  variations 
resulting  from  fluid  and  solid  Earth  tides)  may  not  be  represented  in  the  orbit  propaga¬ 
tion  model  because  inclusion  of  the  full  force  model  is  computationally  expensive.  These 
process  noise  effects,  representing  errors  in  the  motion  models  used  for  orbit  propagation, 
will  impact  the  total  error  in  the  predicted  or  retrodicted  state  of  an  object. 

The  dynamics  of  an  orbiting  object  can  be  represented  in  a  cartesian  Earth-Centered 
Inertial  (ECI)  reference  frame  as 


(  x\ 

(  vx  \ 

°  \ 

y 

vv 

.  0 

z 

Vz 

0 

vx 

ax(x,y,z,vx,vv,vz) 

+ 

Vy 

ay(x,  y,z,vx,vy,vz) 

^y(t) 

W 

\a2(x,y,z,vx,vv,vz)J 

where  px(t),  py(t),  and  pz(t)  correspond  to  errors  in  the  computed  acceleration.  We  as¬ 
sume  that  the  covariance  of  p(t)  is  known  from  estimating  the  magnitude  of  the  neglected 
terms  in  the  force  model. 

For  the  purpose  of  initial  orbit  determination,  the  required  fidelity  of  the  model  in 
use  depends  predominantly  on  the  time  diversity  of  the  initial  set  of  collected  measure¬ 
ments.  If  all  of  the  measurements  to  be  used  are  collected  in  a  relatively  short  window 
(such  as  a  single  pass  of  a  low-Earth  orbiting  satellite),  a  lower- fidelity  model  may  be 
more  appropriate  as  fewer  parameters  must  be  estimated  to  produce  a  solution.  For  orbit 
propagation,  we  would  prefer  to  use  a  sufficiently  high-fidelity  model  that  the  dominant 
error  in  the  result  state  is  the  propagation  of  the  initial  state  error,  and  not  accumulated 
model  errors  across  the  propagation  time.  In  practice  this  is  rarely  achievable;  the  space 
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catalog  encompasses  a  very  large  number  of  objects,  and  high-fidelity  integrating  propaga¬ 
tors  are  too  computationally  expensive  for  most  applications.  Historically,  a  compromise 
has  been  reached  between  extremely  low-fidelity  models  (Keplerian  spherical-Earth  or 
analytic  J2  propagators)  and  high-fidelity  integrators  through  the  use  of  the  Simplified 
General  Perturbations  Satellite  Orbit  Model  4  (SGP4). 

In  this  study,  however,  we  note  that  uncertainty  generation  and  propagation  rely 
somewhat  less  on  the  specifics  of  the  model  in  use.  To  first  order,  the  errors  in  the  initial 
state  estimate  of  an  orbit  are  dominated  by  the  measurement  errors  of  the  radar  sensors  we 
have  considered.  We  are  thus  less  interested  in  the  model  details  and  more  in  the  correct 
representation  of  the  mismatch  between  the  simulator  dynamics  and  the  tracking  model 
dynamics.  For  ease  of  prototyping,  then,  we  have  used  the  simplest  possible  dynamics 
formulation:  that  of  two-body  Keplerian  dynamics  on  a  spherical  Earth.  This  dynamics 
model  is  briefly  described  below: 


F 

1  grav  —  ^  o  •  > 


(2) 


where  G  is  the  universal  gravitational  constant,  m\  is  the  mass  of  the  first  object,  m2  is 
the  mass  of  the  second  object,  r  is  the  difference  in  their  positions,  and  r  is  the  magnitude 
of  that  difference. 

In  the  system  of  interest,  this  causes  an  acceleration  of  the  small  orbiting  body  equal 
to 


Mi 


=  G— rr’ 


xgrav 


(3) 


where  Me  is  the  mass  of  the  Earth,  r  is  the  difference  in  position  between  the  center  of 
the  Earth  and  the  object,  and  r  is  the  distance  between  the  two. 

Note  that  none  of  the  fundamental  methods  to  be  discussed  rely  on  this  model;  the  ap¬ 
proach  applies  equally  well  to  an  arbitrary  model  selection.  The  only  requirement  for  the 
use  of  an  arbitrary  model  is  that  the  process  noise  describing  the  difference  between  the 
model  and  the  expected  real  object  dynamics  be  sufficiently  well  characterized.  We  antic¬ 
ipate  that,  with  some  extensions,  the  SGP4  model  could  be  augmented  with  uncertainty 
propagation  at  an  acceptable  computational  cost. 


3  Approach 

In  order  to  predict  the  target’s  location  at  some  future  point  in  time,  we  must  first  estimate 
its  state  at  the  times  it  is  measured,  and  then  determine  its  most  likely  path  across  the 
intervening  time  interval. 

3.1  Initial  Orbit  Determination 

The  problem  of  determining  the  initial  orbit  is  called  “track  initiation”  in  multiple  target 
tracking  language.  Given  a  sequence  of  measurements  {zi,...,zn}  taken  at 
respectively,  the  goal  is  to  determine  the  position  and  velocity  (or  orbital  elements)  of  the 
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orbiting  target  at  some  time,  usually  within  the  time  spanned  by  the  measurement  times. 
Our  goal  is  to  solve  the  constrained  nonlinear  least  squares  optimization  problem 

n 

Minimize  || zk  -  h  ( x(tk ;  y,  t0))||fl-i  (4) 

k=  1 

o?x 

Subject  to  —  =  /(x,  t) 

at 

*{t0)  =  y. 

Here,  we  are  optimizing  over  the  vector  y  given  a  time  t0 ,  dyanmics  =  f(x,t)  and 
the  above  measurement  sequence.  The  dynamics  may  include  a  general  gravity  model  and 
a  host  of  other  physical  forces  on  the  object  such  as  drag  and  solar  radiation  pressure.  In 
addition,  one  can  include  process  noise  in  the  above  formulation  as  found  in  the  books 
Jazwinski[4]  and  Sage  and  Melsa[9],  or  the  paper  by  Poore  [8]. 

There  are  a  variety  of  optimization  methods  for  solving  the  above  problem,  including 
Gauss-Newton,  full  Newton,  or  hybrid  methods  with  a  globalization  procedure  such  as  a 
line  search  or  trust  region  method.  (The  Levenberg-Marquardt  method  may  be  viewed 
as  a  Gauss-Newton  method  with  a  trust  region  globalization.)  To  initialize  the  nonlinear 
least  squares,  one  could  filter  and  smooth  the  data  or  could  use  three  measurements  to 
approximate  the  position  and  velocity.  The  covariance  of  the  estimate  y  is  identified  in 
the  course  of  solving  the  nonlinear  least  squares  problem  ([9]). 

3.2  Consider  Analysis 

Newton’s  method  yields  a  covariance  matrix  based  on  the  final  solution  and  the  statisti¬ 
cally  formulated  objective  function  (4).  As  a  result  of  both  the  nonlinear  transformation 
between  a  radar  measurement  and  a  target  position,  and  also  the  nonlinear  dynamics  of 
the  system,  the  error  in  the  estimate  does  not  actually  follow  a  Gaussian  distribution  even 
if  the  measurement  errors  are  Gaussian.  In  Cartesian  space,  a  covariance  matrix  can  only 
represent  an  ellipsoidal  uncertainty  region.  The  curved  pancake  from  figure  1  is  flattened 
into  a  highly  eccentric  ellipsoidal  region.  This  region  does  not  accurately  reflect  the  true 
system  uncertainty.  The  two  uncertainty  regions  axe  compared  in  figure  2. 

One  method  of  compensating  for  this  approximation  in  the  covariance  computation 
is  called  a  consider  analysis.  In  a  consider  analysis,  the  covariances  of  the  measurements 
are  increased  to  also  account  for  the  error  caused  by  the  linearization  of  the  nonlinear 
transformation.  Here  is  a  brief  presentation.  We  assume  the  measurement  equation  to 
have  the  form  z  =  h(x,  b)  +  u  where  nu  €  A/^O,  R )  is  the  measurement  noise.  Assuming 
one  has  obtained  an  estimate  of  the  biases  b,  say  b,  one  can  expand  the  measurement 
equation  via  h[x,  b)  «  h(x,  b)  +  Dbh{x ,  6)(e)  where  e  =  b  —  6.  The  modified  measurement 
equation  then  is 

z  =  h(x,  b)  +  Dbh(x,  b)(e)  +  v 

Again,  assuming  e  and  v  to  be  independent  and  t  has  covariance  C,  the  measurement 
error  Rk  in  equation  (4)  is  replaced  by  Rk  +  Dbh(x,  b)CDlh(x ,  b).  The  second  term  gives 
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Figure  2:  The  flattened  uncertainty  region  caused  by  linearizing  the  transformation.  The 
crosses  are  drawn  from  a  linearized  distribution  with  the  same  mean.  The  uncertainty 
region  for  this  distribution  does  not  capture  the  possible  error  of  the  measurements. 

us  the  opportunity  to  account  for  the  residual  biases  and  for  the  accuracy  achievable  in 
the  sometimes  highly  ill-conditioned  transformation  h(x). 

Said  in  another  way,  this  procedure  allows  the  covariance  of  the  final  solution  to 
capture  the  truth,  even  in  cases  where  the  errors  were  actually  larger  than  expected  due 
to  the  transformation  approximation.  This  method  thickens  the  flat  pancake  so  that  the 
ellipsoidal  uncertainty  region  captures  the  actual  uncertainty.  The  effect  is  illustrated  by 
the  much  thicker  ellipsoidal  region  shown  in  figure  3. 

The  Schmidt-Kalman  filter  goes  a  step  beyond  the  above  consider  analysis  and  builds 
the  statistical  cross  correlation  between  the  biases  b  and  the  state  x  and  has  been  shown 
to  be  very  effective  for  improving  covariance  consistency  [7]. 

3.3  Propagating  Uncertainty 

The  equations  of  motion  (1)  can  be  integrated  (along  with  a  covariance  matrix  representing 
the  uncertainty  of  the  state)  using  a  Riccati  equation  to  predict  the  position  of  the  target 
at  a  future  time.  Direct  integration  of  the  equations  of  motion  is  used  when  a  complete 
model  of  the  orbital  dynamics  is  desired,  as  it  allows  for  a  general  acceleration  estimate 
and  error  description  to  be  inserted  into  the  dynamics. 

The  Riccati  equation  for  this  problem  linearizes  the  dynamics  at  each  point  in  the  in¬ 
tegration.  Linear  propagation  provides  a  very  simple  method  of  updating  the  covariances. 
The  dynamics  for  the  covariance  become 

HP 

—  =  F(x)P  +  PF(x)T  +  Q(t), 
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3  APPROACH 


Figure  3:  The  uncertainty  region  has  been  thickened  to  cover  the  actual  errors,  by  using 
consider  analysis.  The  uncertainty  region  represented  by  the  plus  symbols  now  contains 
the  true  errors. 

where  F(x)  is  the  Jacobian  of  the  dynamics  at  the  current  state  x,  and  Q(t)  is  the 
covariance  associated  with  the  random  errors  /xx,  fiy,  and  nz. 

An  alternative  to  integration  is  to  transform  the  cartesian  position  and  velocity  es¬ 
timate  into  orbit  elements.  Orbit  elements  describe  an  elliptical  orbit  and  the  target’s 
position  along  that  orbit;  variations  on  the  classical  Kepler ian  orbit  elements  can  be  used 
to  describe  a  range  of  analytic  solutions  to  the  orbital  dynamics  problem.  The  simplest 
form  of  these  is  that  used  for  a  spherical  Earth  force  model.  These  analytic  methods  are 
much  more  computationally  efficient  than  an  integrating  propagator,  at  the  expense  of 
only  being  able  to  consider  a  subset  of  the  forces  on  a  space  object.  The  transformation 
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3  APPROACH 


from  position,  r,  and  velocity,  v,  into  orbit  elements  is 


IMI.  «  =  IMI, 

r  x  v,  h  =  ||h||, 
k  x  h,  n  =  ||n||, 

^  ((v2_  r)r_^'V^V)’  e  =  Hell« 

hi 

1  -e2’ 

i  <  180° 
h 

— ,  if  ns  >  0  then  Q  <  180° 
n 

-,  if  >  0  then  u>  <  180° 
ne 

if  r  •  v  >  0  then  v  <  180°, 
er 

where  the  orbit  elements  are  a  (semi-major  axis  of  elliptical  orbit),  e  (eccentricity  of 
elliptical  orbit),  i  (inclination  angle  of  the  orbit  from  the  equatorial  plane),  Q  (angle 
of  ascension  through  equatorial  plane  from  the  equinox),  u  (angle  between  perigee  and 
equatorial  passage  within  the  orbital  plane),  and  v  (angle  from  perigee  to  current  target 
position).  For  simple  point  mass  models  of  gravity,  only  the  current  position  along  the 
orbit  v  needs  to  be  updated.  The  relationship  between  the  time  that  has  passed  since  the 
last  time  the  target  passed  perigee  and  the  angle  v  is  described  by  Kepler’s  equation  [1]. 


r  = 

h  = 
n  = 

e  = 

a  = 

cos(i)  = 
cos(fi)  = 
cos  (a;)  = 
cos(i/)  = 


cos  (E)  = 


e  cos(i/) 

1  +  e  cos(u) 


A  t  = 


1  —  {E  -  esin(E)) 


This  second  method  of  propagation  does  not  inherently  allow  for  appropriate  propa¬ 
gation  of  uncertainty.  However,  because  the  propagation  is  extremely  quick,  a  particle 
transform  can  be  employed.  In  a  particle  transform,  many  pseudo-random  values  from 
an  initial  distribution  are  drawn,  and  each  one  is  propagated  the  appropriate  amount  of 
time.  If  the  dynamics  model  has  random  errors,  small  appropriate  perturbations  can  be 
added  to  each  particle  during  the  propagation.  After  propagation  the  values  represent 
the  distribution  at  the  later  time.  Figure  4  illustrates  the  results  of  the  method. 
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4  SIMULATION  RESULTS 


Figure  4:  The  thick  line  represents  the  mean  orbit  of  the  particles  (black  dots).  The  true 
orbit  illustrated  by  the  thin  blue  line  lies  within  the  cloud  of  particles.  While  there  is 
error,  the  covariance  of  the  particles  in  the  cloud  clearly  captures  the  blue  star  where  the 
target  is  truly  located. 

4  Simulation  Results 

4.1  Normalized  Chi-Squared  Metric 

In  order  to  test  whether  the  covariance  matrices  developed  using  these  approaches  were 
appropriate,  we  examine  covariance  consistency.  The  mean  squared  error  of  an  unbiased 
estimator  that  follows  a  Gaussian  distribution  with  unit  variance  follows  a  chi-squared 
distribution  with  one  degree  of  freedom.  This  distribution  has  an  expected  value  of  one. 
Furthermore,  the  mean  squared  error  of  an  estimate  divided  by  the  variance  of  the  estimate 
also  follows  this  distribution.  Analogously,  the  quantity 

(x  -  x*)TP~\x  -  x*),  (5) 

where  x  is  drawn  from  an  m-dimensional  Gaussian  distribution  with  mean  x*  and  co- 
variance  P~l  will  have  a  chi-squared  distribution  with  m  degrees  of  freedom,  which  has 
a  mean  of  m.  We  normalize  by  dividing  out  the  dimension  m  and  the  mean  once  again 
equals  one. 

In  order  to  determine  if  the  covariance  matrices  developed  using  the  methods  described 
is  an  adequate  representation  of  the  errors  in  the  estimates,  we  can  examine  the  normalized 
chi-squared  metric.  If  this  value  is  far  greater  than  one  on  average,  then  the  covariance 
matrix  is  actually  to  small  to  capture  the  errors.  If  the  value  is  less  than  one  on  average, 
then  the  covariance  matrix  is  too  large  on  average,  and  the  actual  errors  are  smaller  than 
it  indicates. 
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4.2  Scenario  Description 

A  single  target  in  orbit  is  produced  following  a  Keplerian  trajectory.  It  passes  through  the 
field  of  view  of  a  radar  platform  that  is  sitting  on  the  Earth.  The  radar  generates  mea¬ 
surements  of  a  target  by  first  transforming  the  true  target  position  into  its  measurement 
space,  and  then  adding  independent  Gaussian  pseudo-random  errors  to  these  quantities. 
The  covariances  of  these  quantities  are  recorded  along  with  the  measurements. 

We  run  the  simulation  in  two  different  modes.  In  the  first  mode,  no  consider  analysis 
is  performed  on  the  measurements;  the  measurements  are  passed  directly  into  the  batch 
estimator.  In  the  second  mode,  a  linearization  compensation  error  covariance  of  50  m2  in 
each  cartesian  coordinate  is  added  to  the  measurement  covariance  through  the  appropriate 
transformation.  This  additional  uncertainty  keeps  the  uncertainty  region  described  by  the 
covariance  from  being  thin  and  flat,  when  it  should  be  curved. 

The  nonlinear  least  squares  problem  (4)  is  solved  to  produce  an  estimate,  x,  and  a 
covariance  matrix  Px.  The  covariance  consistency  metric  is  computed  for  this  estimate. 
The  next  step  is  to  propagate  the  track  states  and  their  uncertainties  for  12  hours.  To 
do  this  we  use  the  two  different  methods  described:  integrate  the  equations  of  motion 
using  the  Ricatti  equations  to  carry  the  covariance,  or  use  a  particle  transform  into  orbit 
elements,  and  perform  the  propagation  in  that  space.  In  these  simulations  the  true  target 
trajectory  was  computed  using  a  spherical  Earth  model,  without  any  additional  forces. 
Therefore,  the  process  noise  term  Q(t )  was  zero,  and  no  random  changes  to  the  orbit 
element  particles  were  added  during  propagation.  The  particle  transformation  used  a 
sample  size  of  1000.  After  the  propagation,  the  covariance  consistency  metric  is  computed, 
using  an  integrated  spherical  Earth  motion  model  to  produce  the  true  state  at  the  later 
time. 


4.3  Data 

The  number  of  measurements,  and  the  amount  of  time  they  span  greatly  affects  the 
ability  of  the  batch  estimation  process  to  correctly  estimate  the  target  state.  With  few 
measurements  or  a  very  short  span  of  time,  the  estimator  produces  poor  estimates.  With 
many  measurements  spread  over  a  long  period  of  time,  the  estimator  provides  useful 
estimates.  The  driving  force  behind  this  phenomenon  is  the  velocity  portion  of  the  original 
estimate.  When  only  a  few  measurements  are  available,  and  they’re  all  at  nearly  the  same 
time,  the  measurement  errors  prevent  an  accurate  estimate  of  the  velocity  from  being 
obtained.  If  the  measurements  are  widely  spaced,  or  enough  measurements  are  taken, 
this  effect  is  mitigated. 

The  scenario  described  was  run  with  between  5  and  40  measurements  with  average 
total  span  of  between  4  and  195  seconds.  Several  variations  of  number  of  measurements 
and  time  spacing  of  measurements  were  simulated.  Each  variation  was  performed  in  500 
Monte  Carlo  runs  in  each  mode  (with  or  without  consider  analysis).  The  covariance 
consistency  metric  was  averaged  over  these  500  runs. 

After  performing  these  experiments,  we  notice  that  the  integration  of  the  covariance 
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4  SIMULATION  RESULTS 


Table  1:  Mean  Normalized  Chi-Squared  Statistic  from  500  MC  Runs  with  Consider  Anal¬ 
ysis 


Measurements 

Spacing 

Span 

To  ECI  x2 

Integration  x2 

T0  OE  x2 

Particle  x2 

5 

1 

4 

2.5060 

- 

1.4423 

1.3429 

5 

5 

20 

2.4973 

- 

1.1934 

1.0300 

5 

10 

40 

2.4037 

- 

1.0513 

1.0545 

5 

20 

80 

2.5733 

- 

1.1189 

1.1941 

5 

40 

160 

1.7023 

- 

.70684 

.79285 

10 

1 

9 

1.6321 

- 

1.3547 

.93299 

10 

5 

45 

1.6691 

- 

.98633 

.96807 

10 

10 

90 

1.5266 

- 

.92922 

.98162 

10 

20 

180 

1.3378 

- 

.79309 

.86993 

20 

1 

19 

1.2581 

- 

.97881 

.92262 

20 

5 

95 

1.1174 

- 

.84290 

.88579 

20 

10 

190 

1.0001 

- 

.76438 

.86928 

40 

1 

39 

.96285 

- 

.85410 

.80293 

40 

5 

195 

.82877 

- 

.70033 

.76643 

Table  2:  Mean  Normalized  Chi-Squared  Statistic  from  500  MC  Runs  without  Consider 
Analysis 


Measurements 

Spacing 

Span 

%  ECI  x2 

Integration  x2 

%  OE  x2 

Particle  x2 

5 

1 

4 

4.3587 

- 

2.3315 

2.1651 

5 

5 

20 

3.8840 

- 

1.4460 

1.2805 

5 

10 

40 

3.0420 

- 

1.1065 

1.1001 

5 

20 

80 

3.7265 

- 

1.1288 

1.1843 

5 

40 

160 

2.8620 

- 

1.0555 

1.0787 

10 

1 

9 

2.4654 

- 

1.7041 

1.1698 

10 

5 

45 

2.2713 

- 

1.1022 

1.1202 

10 

10 

90 

2.4220 

- 

1.2131 

1.2611 

10 

20 

180 

2.0620 

- 

1.0976 

1.1622 

20 

1 

19 

1.7244 

- 

1.0809 

1.0095 

20 

5 

95 

1.5003 

- 

.99134 

1.0131 

20 

10 

190 

1.3443 

- 

.98115 

.97296 

40 

1 

39 

1.4096 

- 

.99437 

1.0459 

40 

5 

195 

1.2076 

- 

1.0011 

1.0612 
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using  the  Riccati  equation  cannot  be  performed  in  a  stable  manner  with  zero  process 
noise.  The  integration  quickly  deteriorates  the  condition  number  of  the  matrix,  until  it 
is  numerically  singular,  and  no  longer  has  a  valuable  meaning.  Without  the  covariance 
matrix,  the  chi-squared  consistency  metric  cannot  be  computed. 

Based  on  tables  1  and  2,  we  can  see  that  if  enough  measurements  are  available  spread 
over  a  long  enough  time  period,  then  the  estimates  represented  as  orbital  elements  can 
achieve  consistent  covariances,  without  the  use  of  consider  analysis.  This  is  not  the  case 
for  the  ECI  estimates,  however.  The  elliptical  uncertainty  region  represented  by  the  ECI 
covariance  matrix  is  the  wrong  shape.  If  it  is  to  capture  the  true  errors  in  a  consistent 
manner,  a  consider  analysis  performed  for  the  specific  number  of  measurements,  timing, 
and  geometry  must  be  performed.  The  additional  uncertainty  of  50  m2  in  each  direc¬ 
tion  did  make  the  scenario  with  20  measurements  taken  one  every  10  seconds  consistent. 
However,  this  consider  value  was  excessively  large  for  the  scenario  with  40  widely  spread 
measurements,  and  insufficient  for  the  scenarios  with  only  10  measurements. 

5  Conclusion 

The  goals  of  this  project  were  twofold:  (1)  given  a  sequence  of  radar  measurements, 
establish  an  object’s  position  and  velocity  at  a  specified  time  and  develop  a  consistent 
covariance  matrix  that  represents  the  uncertainty  in  the  initial  state,  and  (2)  determine  a 
method  of  propagating  this  state  and  its  uncertainty  that  maintains  a  consistent  uncer¬ 
tainty  measure  for  a  period  of  hours  or  days. 

We  achieved  the  first  of  these  goals  by  using  Newton’s  method  combined  with  a  con¬ 
sider  analysis  as  demonstrated  in  tables  1  and  2.  One  should  the  bias  covariances  arising 
from  a  bias  estimator  to  better  judge  the  effectiveness  of  the  algorithm. 

The  second  goal  was  attained  by  using  a  particle  transform  to  describe  the  states  in 
orbital  element  space.  This  is  demonstrated  in  figure  4,  as  well  as  with  the  metrics  in 
tables  1  and  2.  We  also  attempted  to  use  the  Riccati  equation  in  ECI,  but  this  led  to 
an  unrealistic  covariance  procedure.  Perhaps  the  use  of  the  Riccati  equation  and  then 
conversion  back  to  orbital  element  space  may  be  more  appropriate. 

In  future  work,  a  Schmidt-Kalman  filter-smoother  might  be  used  instead  of  Newton’s 
method  to  appropriately  take  into  account  the  statistical  cross  correlation  between  biases 
and  the  state  to  be  estimated.  The  inclusion  of  higher  order  gravity  model  and  such  forces 
as  atmospheric  drag  and  solar  radiation  pressure  ([6])  in  the  dynamics.  The  appropriate 
distribution  of  perturbations  to  be  added  to  particles  would  need  to  be  researched.  Fi¬ 
nally,  we  would  like  to  test  these  methods  using  different  sensor  models,  such  as  optical 
telescopes. 
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